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We study the behaviour of double-stranded RNA under twist and tension using oxRNA, a re¬ 
cently developed coarse-grained model of RNA. Introducing explicit salt-dependence into the 
model allows us to directly compare our results to data from recent single-molecule experiments. 
The model reproduces extension curves as a function of twist and stretching force, including 
the buckling transition and the behaviour of plectoneme structures. For negative supercoiling, 
we predict denaturation bubble formation in plectoneme end-loops, suggesting preferential plec¬ 
toneme localisation in weak base sequences. OxRNA exhibits a positive twist-stretch coupling 
constant, in agreement with recent experimental observations. 


I. INTRODUCTION 

Due to their importance in the storage and pro¬ 
cessing of genetic information, nucleic acids play a 
fundamental role in many biological processe s su ch 
as transcription, translation and replication.^!^ In 
their double stranded (ds) form, DNA and RNA 
adopt a helical geometry. While dsDNA typi¬ 
cally forms a B-helix, dsRNA adopts an A-helical 
form, which is wider, has a smaller pitch and 
bases that are inclined with respect to the heli¬ 
cal axis.l^ Double-stranded DNA and RNA exhibit 
complex mechanical behaviour that is important 
in many biomechanical conte xts, such as genome 
organisatioi^ virus packagin^^ and nucleosome 

positioning.13 

Moreover, both DNA^^ and more recently 
RNA^i^ have emerged as versatile building materi¬ 
als on the nanoscale. Driven by these wide-ranging 
applications, the mechanical properties of nucleic 
acids have been studied with increasing precision 
on a single-molecule level.l^ While the mechanical 
behaviour of dsDNA has been widely characterised 
using molecular tweez er as says dsRNA has 
received less attention.l^SIIIl The first comprehen¬ 
sive experimental study of the twisting and stretch¬ 
ing behaviour of dsRNA was only recently carried 
out by Lipfert and co-workers.!^ 

Correspondin gly, th eoretical work usi ng at om- 
istic simulationscontinuum model^^^^^fil 
coarse-grained simulationd^IHSai j^gg centered on 
modelling the properties of torsionally stressed 
DNA. There have been far fewer studies of su¬ 
percoiled dsRNA, although theoretical investiga¬ 
tions exist using atomistic simulation^^ni g^gt the 
HelixMC package, which uses a base-pair-level de¬ 
scription of the molecule.!^ 


Here, we study the behaviour of super- 
coiled dsRNA using a salt-dependent extension 
of oxRNA, a recen tly developed nucleotide-level 
model of RNA.f^^^The model is developed to cap¬ 
ture the structural, mechanical and thermodynam¬ 
ical properties of both single-stranded and double- 
stranded RNA and was previously used to study 
RNA hairpin unzipping, the thermodynamics of 
pseudoknot folding, kissing complex for mation and 
toehold-mediated strand displacement .121133] 
coarse-graining methodology of oxRNA allows us 
to capture the effects of double-strand denatura¬ 
tion, which are not accessible in continuum or 
basepair-level models. Likewise, the computa¬ 
tional efficiency gained by the coarse-graining al¬ 
lows us to access time scales and system sizes 
relevant to the physics of double-strand buck¬ 
ling and denaturation, which are currently beyond 
the scope of all-atom molecular dynamics simula¬ 
tions. We previ ously used a coarse-grained model 
of DNA, oxDNA! 22E11^ to study the supercoiling of 
dsDNA and o btaine d good agreement with experi¬ 
mental results.!2S125lin this work, we use oxRNA to 
directly compare to a recent experimental study of 
dsRNA supercoiling.!^ 

This paper is organised as follows. First, we 
briefly describe an extension of the oxRNA model 
to include a salt-dependent parameterisation. We 
then compare the model prediction to recent mea¬ 
surements of the end-to-end distance and torque 
response of dsRNA as a function of imposed 
stretching force and superhelical density.!^ We ex¬ 
tract parameters characterising the twisting, bend¬ 
ing and extensional behaviour of the molecule. The 
results of our simulations are in reasonable agree¬ 
ment with experimental data. For negative super¬ 
coiling and intermediate stretching forces, we ob- 
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FIG. 1. A schematic representation of (a) an A-RNA 
helix as represented by the oxRNA model and (b) 
the attractive interactions in oxRNA. The lines in (b) 
schematically show the interactions between the nu¬ 
cleotides: Hydrogen bonding (Vh.b.), stacking (I4tack), 
cross-stacking (Fcroas st.) between a nucleotide and the 
nucleotide that is the 3' neighbour of the directly oppo¬ 
site nucleotide and coaxial stacking (I4oaxiai at.)- The 
nucleotides also interact with excluded-volume inter¬ 
actions and electrostatic interactions, which are not 
shown. 


serve denaturation bubble formation localised in 
plectoneme end-loops, similarly to what was found 
in a previous work on DNA plectonemes using a 
related modelling approach for DNA.I^ 


II. OXRNA MODEL WITH SALT-DEPENDENT 
INTERACTION 


OxRNA represents each nucleotide as a single 
rigid body with multiple interaction sites. The 
rigid bodies interact with effective anisotropic in¬ 
teractions that are designed to capture the overall 
thermodynamic and structural consequences of the 
base-pairing, stacking and backbone interactions, 
as schematically shown in Fig. The potential of 
the oxRNA model is 


hoxRNA — ^ ^ ^Rbackbone “b l^tack T 

(d> 

+ + Rcross st. T l^XC 

“b fAoaxial st. “b V^lectrostatic) j (1) 


where the first sum runs over all pairs of nu¬ 
cleotides which are nearest neighbours on the same 
strand and the second sum runs over all other 
pairs. A detailed description of the interactions 
and their parameterisation is provided in Ref. 1311 
with the exception of I4iectrostatic which is newly 
introduced to explicitly capture salt-dependent ef¬ 
fects. This term is isotropic and is centred on the 
backbone site of each nucleotide. The functional 
form of the potential is based on Debye-Hiickel the¬ 
ory, where we further introduce a cutoff at a finite 


distance. We use the Debye-Hiickel length for wa¬ 
ter and treat the strength of the effective negative 
charge on the backbone site as a parameter, which 
we fit to reproduce the melting temperatures of 
duplexes of lengths 5, 6, 7, 8, 10 and 12 at salt con¬ 
centrations varying from 0.1 M to IM. To obtain 
the melting temperatures to which we fit, we use 
the averaged nearest-neighbour model of Turner et 
extended with a salt-dependent free-energy 
correction inferred from hairpin unzipping experi¬ 
ments at varying salt conditions.^ We employ the 
htting procedure based on thermodynamic integra¬ 
tion, as detailed in Ref.|3Hl We provide further de¬ 
tails of the functional form of Wiectrostatic and its 
parameterisation in the Supplementary Material.l^ 

The backbone interaction, Vbackbone) is an 
isotropic FENE spring potential that is used to 
mimic the covalent bonds in the RNA backbone 
that constrain the intramolecular distance be¬ 
tween neighbouring nucleotides. The nucleotides 
further have repulsive excluded-volume interac¬ 
tions 14xc and 14 xc ftiat depend on the dis¬ 
tance between their interaction sites, namely 
the backbone-backbone, stacking-stacking and 
stacking-backbone distances. The excluded- 
volume interactions ensure that strands cannot 
overlap, or pass through each other in a dynam¬ 
ical simulation. 

The duplex is stabilised by hydrogen bond¬ 
ing (Rh.b.), stacking (14tack) and cross-stacking 
(Rcross st.) interactions. These potentials are 
anisotropic and depend on the distance between 
the relevant interaction sites as well as the mu¬ 
tual orientations of the nucleotides. The hydrogen¬ 
bonding term Vh.b. captures the stabilising in¬ 
teractions between complementary Watson-Crick 
(AU and GC) and wobble (GU) base pairs, while 
Rstack mimics the favourable interaction between 
adjacent bases on the same strand. The strength 
of Vh.b. and 14tack is sequence-dependent, i.e. de¬ 
pends on the identity of the interacting bases. 

The cross-stacking potential, Wross st., is de¬ 
signed to capture the interactions between diag¬ 
onally opposite bases in a duplex and has its min¬ 
imum when the distance and mutual orientation 
between nucleotides corresponds to that for a nu¬ 
cleotide and the 3' neighbour of the directly oppo¬ 
site nucleotide in an A-form helix. This interaction 
has been parameterised to capture the stabilisation 
of an RNA duplex by a 3' overhang. 

The coaxial stacking potential Rcoaxiai st. repre¬ 
sents the stacking interaction between nucleotides 
that are not nearest neighbours on the same 
strand. 

In this work, we use the average-base param¬ 
eterisation of oxRNA, which only allows for spe- 
cihc formation of AU and GG Watson-Crick base 
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pairs. Hydrogen-bonding energies between com¬ 
plementary base-pairs and stacking energies are set 
to identical, average strengths. This choice allows 
us to focus on the generic properties of RNA dou¬ 
ble strands, which are independent of specific se¬ 
quence properties. Parameters are fitted to repro¬ 
duce the thermodynamics of hairpins and duplexes 
averaged over all possible combinations of Watson- 
Crick base pair steps, as predicted by the model 
of Turner and collaborators.!^ We note that the 
model cannot reproduce tertiary structure contacts 
such as ribose zippers or Hoogsteen base pairs, but 
we do not anticipate that these non-canonical in¬ 
teractions will be relevant for the modelling of the 
behaviour observed in Ref. [221 


III. SIMULATION METHODS 

The results reported in this work were obtained 
from molecular dynamics simulations of oxRNA 
using an Andersen-like thermostat (described in 
the appendix of Ref. 001) at 300 K using both the 
CPU and GPU implementation of the model.!^ We 
intentionally set the diffusion constant artificially 
high to speed-up convergence of the simulations 
to equilibrium. In particular, we used the trans¬ 
lational diffusion constant D = 5.8 x 10“^m^s“^, 
which corresponds to a diffusion constant of 2.1 x 
10“®m^s“^ for a 14-mer and is about two orders 
of magnitude more than the experimentally mea¬ 
sured I?exp = 0.92x 10“^° m^s“^.!^The simulation 
time step was set to 1.22 x 10“^"^ s. 

We simulated 600-bp dsRNA molecules using an 
average-base parameterisation of oxRNA that in¬ 
cludes base-pair specificity, but ignores sequence- 
dependent variations in interaction energies The 
duplex was set-up as a homogenously twisted he¬ 
lix with a desired superhelical density and pre¬ 
equilibrated for a simulation time of at least 1 fis. 
Simulations were then run for at least 8 /iS of sim¬ 
ulation time. The superhelical density is defined 
as (7 = pq/p — 1, where p is the imposed pitch 
and pq is the equilibrium pitch of dsRNA when 
no stress is applied. To keep superhelical densi¬ 
ties constant during a simulation run, strand ends 
were fixed in two-dimensional harmonic traps and 
the strands prevented from passing around their 
own ends, as described in detail in the Supplemen¬ 
tary Material.!^ The resulting setup of the dsRNA 
systems subject to linear and torsional stress is il¬ 
lustrated schematically in Fig. 

To match the experimental conditions of Ref. [221 
all simulations were run at a monovalent salt con¬ 
centration of 100 mM. 



FIG. 2. Schematic of the simulation setup used in this 
work. A-helical dsRNA strands are subjected to tor¬ 
sional stress by fixing a constant superhelical density 
a and exerting a stretching force F to the strand ends. 
Repulsive planes tagged to the strand ends (indicated 
in grey) ensure the superhelical density remains con¬ 
stant by preventing the duplex from passing around its 
own ends. The conhguration shown was obtained in a 
simulation a,t a — -1-0.08 and F = 3.0 pN. Under these 
conditions, a plectoneme forms leading to significant 
shortening of the end-to-end extension. 


IV. RESULTS 


Superhelical stress can be stored in dsDNA and 
dsRNA by both twisting and writhing. For small 
values of supercoiling, the torsional energy of the 
system grows until a buckling superhelical density 
(Tb is reached, at which it becomes more favourable 
for the system to form writhed structures known 
as plectonemes (see Fig.[^ where the supercoiling 
energy is stored in bending rather than twisting.!^ 
Writhing, which results in a shortening of the 
molecule end-to-end distance, is disfavoured by ap¬ 
plying an external stretching force. As described in 
more detail below, our model exhibits this generic 
behaviour, as expected for a twist-storing polymer 
with finite bending persistence length. Here we 
compare the behaviour of our model to experimen¬ 
tal data of Lipfert and co-workers.!^ 
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FIG. 3. Mechanical behaviour of 600-bp dsRNA under torsional stress in the oxRNA model at lOOmM monovalent 
salt, together with experimental data obtained for a 4.2 kbp-system from Ref. 1221 (a) “Hat curves” showing 
the end-to-end extension of the system. Arrows indicate the values of a for which buckling is expected from 
Eq. m Error bars (standard deviations) indicate the magnitude of thermal fluctuations rather than measurement 
uncertainties, (b) Torque response of the dsRNA strand, showing a linear regime at low twist, followed by a 
constant-torque regime after buckling, (c) Postbuckling slopes measured from the hat curves in (a) and a fit 
of simulation results to the analytical model of Ref. 1441 leading to a torsional stiffness of the plectonemic state 
ToxRna = 22 nm. (d) Fit of the torsional stiffness measured in simulations to a Moroz-Nelson model.^ 


A. Force-extension response at varying superhelical 
densities 

We first study the end-to-end extension of a 600- 
bp dsRNA as a function of superhelical density a 
and stretching force F. For a given dsRNA with 
imposed a and F, we run a molecular dynamics 
simulation, as described in Section |III[ and mea¬ 
sure the end-to-end distance between the first and 
the last base pairs of the duplex. The results are 
shown in Fig. [^a). When the superhelical density 
of the dsRNA molecule in our model is increased, 
its end-to-end extension initially changes little, un¬ 
til a buckling point is reached at which it is ther¬ 
modynamically more favourable for the system to 
bend into a plectonemic structure than to further 
twist. For stretching forces F > 2pN, the exten¬ 
sion curves become asymmetric, as denaturation 
rather than plectoneme formation occurs for neg¬ 


ative supercoiling (Fig. [^. 

Comparison of our simulation results to the re¬ 
cent experimental data of Ref. (included in 
Fig. I^a)) shows good agreement for the buckling 
superhelical densities, post-buckling slopes, and 
the onset of double-strand melting, indicating that 
the overall behaviour of dsRNA subject to twist 
and stretching force is well reproduced by oxRNA. 

Nevertheless, oxRNA still buckles under positive 
supercoiling for stretching forces above 5 pN, while 
no buckling was observed in experiment above such 
a force.^It was proposed that overwound dsRNA 
above 5 pN changes its conformation to a “P-RNA” 
state that is similar to the P-DNA structure of ds- 
DNA, which is characterised by interwound sugar- 
phosphate backbones with exposed bases.l^Such a 
structure is not observed with oxRNA under these 
conditions. 

Compared to the experimental data, simulated 
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FIG. 4. Strand configurations observed at a stretching 
force F = 2.0 pN for a — +0.09 (left) and cr = —0.09 
(right). At negative superhelical density, a denatura- 
tion bubble of size 10 bp is observed in the end-loop of 
the plectoneme, analogous to the behaviour predicted 
for dsDNA in Ref. 1281 Enlarged structures show the 
microscopic configuration of the plectoneme end-loop, 
where denatured bases are coloured green. 


dsRNA molecules show a larger relative end-to- 
end extension. This is due to the relatively low 
value of the extension modulus AToxRNA ~ 116 pN 
in oxRNA, which is significantly lower than the ex¬ 
perimental value of ATexp ~ 350 pN.I^ However, for 
sufficiently low forces, the buckling behaviour of 
the strand is expected to be only minorly affected 
by this discrepancy. 

As was done in the experimental study, we 
further determined the twist-stretch coupling by 
measuring the slope of the end-to-end exten¬ 
sion curve at low superhelical densities (—0.02 < 
cr < +0.025) and high stretching force F = 
6.0 pN. We obtain a twist-stretch coupling of 
(dAL/dLfc)oxRNA = —0.72 nm/turn, which is 

to be compared to an experimental value of 
{dAL/dLk)e^p = —0.85 nm/turn.l^ Thus, oxRNA 
qualitatively reproduces the positive twist-stretch 
coupling observed for RNA. 

To further quantify the mechanical behaviour 
of oxRNA, we measured the slopes of the exten¬ 
sion curves in the postbuckling regime (shown in 
Fig-iu: c)). We note that the values obtained are 
sensitive to the selection of points included in the 


fit of the postbuckling slope, as indicated by the 
error bars in Fig. [^c). The fitting procedure is de¬ 
scribed in detail in the Supplementary Material.l^ 
Again, approximate agreement with experimen¬ 
tal values is found. When fitting to a thermody¬ 
namic model of the plectonemic phase,qualita¬ 
tively similar but more pronounced systematic de¬ 
viations occur compared to the experimental data, 
as shown in Fig. ic). At least part of the dis¬ 
crepancies may be due to finite size effects in the 
simulated 600-bp system, which approximates the 
thermodynamic limit le ss wel l than the 4.2-kbp ex¬ 
perimental system does.l2S12Sl 

As in a recent study on dsDNA using the oxDNA 
model,we observed localisation of double-strand 
denaturations in the end-loop of plectoneme struc¬ 
tures (see Figs. 1^ andfor cr < 0 and intermediate 
stretching force F « 2pl4^. In this configura¬ 
tion, the enthalpic cost for opening the bubble is 
partially compensated by the lower bending energy 
of a plectoneme end-loop containing a denatura- 
tion bubble; the bubble also reduces the torsional 
stress by absorbing negative twist. 

We note however that the prevalence of these 
bubbles co-localised in the end loops of the plec- 
tonemes is reduced compared to the analogous 
setup in dsDNA. Primarily, this difference may be 
attributed to the stronger base-pairing of Watson- 
Crick base pairs in RNA compared to DNA,I^ 
making bubble opening in stressed parts of the 
strand more enthalpically costly. More subtle ef¬ 
fects, such as differences between the A-form heli¬ 
cal geometry of dsRNA and the B-form helical ge¬ 
ometry of dsDNA, as well as details of the model 
of screened electrostatic interactions may further 
contribute to the differences observed. We also 
note that the simulations presented in this work 
were done at 0.1 M monovalent salt rather than 
the 0.5 M used in Ref. [28] with oxDNA. At vari¬ 
ance with the dsDNA case,l^ we observed double 
strand denaturation only for cr < 0 (Fig. |^, while 
no significant denaturation occurred for positive 
supercoiling at the stretching forces studied in this 
work. This may again be explained by the stronger 
base-pairing free-energy in dsRNA. Force-induced 
melting of the duplex is expected also for cr > 0 
at forces significantly higher than the ones used in 
this work or in experimental assays.l^ 

In this work, we used an average-base parame- 
terisation of oxRNA. However, stable occurrence 
of a tip-bubble plectoneme state for cr < 0 and in¬ 
termediate F « 2 pN suggests that in a sequence- 
dependent scenario, the centres of the plectonemes 
will be primarily localised to AU-rich regions of the 
strand at these conditions, because their weaker 
base paring reduces the cost of bubble formation; 
this mechanism is described in detail for dsDNA 
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in Ref. We note that the occurrence of co¬ 
localised denaturation and writhing is the conse¬ 
quence of the elastic properties of a chiral, semi- 
flexible polymer combined with the possibility for 
the double strand to denature, and is therefore ex¬ 
pected to be a robust phenomenon that is largely 
independent of detailed microscopic properties of 
the molecule. 


B. Torque response and mechanical parameters of 
dsRNA 


We further quantify the properties of dsRNA by 
studying the torque response of molecule at differ¬ 
ent superhelical densities and forces. The torque 
response of the simulated system to imposed super¬ 
helical density is shown in Fig. ib), and compared 
to the corresponding experimental data. Overall, 
we observe fair agreement with the corresponding 
experimental values. 

For small absolute values of the superhelical den¬ 
sity, the torque response of the system grows lin¬ 
early with u. In this regime, the effective torsional 
rigidity of the system corresponds to the slope of 
the torque response curve. The bending and twist 
persistence lengths Aq and Cq can be determined 
by fitting the effective torsional rigidities CcS to a 
model due to Moroz and Nelsorl^(see Fig. [^d)): 



The fits yield Ao,oxRNA = 32 nm and C'o,oxRNA = 
79 nm for the simulated system. Both values are of 
the correct order of magnitude, but lie below the 
values Ao,exp = 57 nm and Co,exp = 100 nm deter¬ 
mined from the experimental systems in Ref. 1221 
The difficulty of correctly reproducing the persis¬ 
tence length in a coarse-grained model of RNA has 
been noted before,!^ and has also affected other 
coarse-grained modelling approaches.!^ However, 
as the relative deviations in the elastic persistence 
lengths are of similar magnitude, we expect prop¬ 
erties that only depend on the ratio of twisting and 
bending energies, such as twist-induced double¬ 
strand buckling to be reproduced more accurately 
by our model than properties that depend on their 
values separately. 

As |cr| is increased, a buckling point is reached 
at which the system forms a plectoneme structure, 
thus absorbing supercoiling by writhing rather 
than further twisting, as discussed previously. 
Buckling occurs once the superhelical density ex¬ 
ceeds a critical value cr{,, which is set by the ratio 
of Co and Aq. The critical superhelical density can 


be estimated bjS^IlIl 


o-fe 


2FAo ropo 
fegT 27rCo’ 


( 3 ) 


where F is the applied stretching force, and cq = 
0.28 nm and po = 11.14 bp are the equilibrium rise 
and pitch of the dsRNA helix, respectively.!^ Using 
the persistence length values obtained by fitting to 
Eq.m tTfe can be predicted from Eq. As indicated 
by arrows in Fig.j^a), the critical superhelical den¬ 
sities obtained in this way are consistent with the 
buckling behaviour observed in simulations. They 
are also consistent with experiment, although it 
should be kept in mind that part of the accuracy 
arises because both Aq and Co are under-estimated 
in oxRNA. Properties which depend on just one of 
these constants will likely agree less well with ex¬ 
periment. 

For low stretching forces, we furthermore ob¬ 
serve a torque “overshoot” (the increase of torque 
before reaching the saturated regime with in¬ 
creased superhelical density) upon buckling, as 
was found experimentally for both DNA^i^ and 
RNA.!^ This overshoot is due to the need to nu¬ 
cleate the end loop of the plectoneme and its mag¬ 
nitude is set by the difference between the free- 
energy cost of forming the plectoneme end-loop 
and the free-energy cost of adding one superhe¬ 
lical turn to an existing plectoneme.^ Decreas¬ 
ing the solvent ionic strength and hence increas¬ 
ing the electrostatic strand repulsion is expected 
to change the free-energy of the relatively large 
end-loop less than that of additional, more tightly 
wound plectoneme turns. Therefore, a reduction of 
the overshoot with decreasing salt concentration 
is expected.!^^ Consistently, we observe a smaller 
overshoot compared to analogous simulations of 
DNA at 500 mM monovalent salt concentration.!^ 

The mechanical parameters of our model ob¬ 
tained so far can be used to derive the torsional 
stiffness of the plectonemic state P by fitting to 
an analytical model introduced by Marko,!^as ex¬ 
plained in detail in the Supplementary MateriaP^ 
(see Fig. [^c)). While trend and order of magni¬ 
tude agree, the simulation results deviate from the 
theoretical prediction due to finite size effects. Ex¬ 
perimental measurements of Ref. from a 4.2-kbp 
dsRNA system show a qualitatively similar devia¬ 
tion from the analytical model, suggesting that at 
higher forces the postbuckling slopes are slightly 
higher than predicted by the analytical model. 

We summarize the mechanical parameters of 
dsRNA inferred in this study for the coarse¬ 
grained model at a monovalent salt concentration 
of 100 mM in Table |Tj along with the correspond¬ 
ing values determined from experiments. In or¬ 
der to be consistent with common experimental 
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FIG. 5. Characterising the tip-bubble state in dsRNA plectonemes: (a) Frequency of plectoneme, co-localised 
tip-bubble and pure bubble states as a function of applied stretching force at cr = —0.09. (b) Representative RNA 
configurations at a = —0.09 and different stretching forces, with bubble positions indicated by green arrows. At 
F = 1.0 pN, no stable strand denaturation occurs, while a 5-bp denaturation bubble localised at the plectoneme 
tip is found at F = 2.0 pN, and a pure, writhed denaturation bubble of size 20 bp occurs at F = 3.0 pN. 


TABLE I. Summary of mechanical parameters deter¬ 
mined for dsRNA at 100-150 mM monovalent salt con¬ 
centration. 


Parameter 


oxRNA experiment 


Bending persistence length Aq [nm] 32 

Torsional persistence length Cq [nm] 79 

Torsional stiffness of plectonemes P [nm] 22 
Extension modulus K [pN] 116 

Equilibrium helical pitch pQ [bp] 11.14 

Equilibrium twist angle 9q [deg] 33.3 

Twist-stretch coupling 

dAL/dLk [nm/turn] -0.72 


57 gpomil] 

iQjm 

2(P1 

350-50dSIil 

10.7- li3 

32.7- 33. P 
-0.8^ 


protocols,^ the equilibrium twist angle Oq and the 
corresponding pitch po = 27r/do were obtained by 
demanding that the overall torque r(F, d) exerted 
on the strand by the traps vanish in a system with 
that twist angle: r(F, 6*0) = 0. 


V. SUMMARY AND CONCLUSIONS 


grained simulations, both for oxRNA,!^ as well as 
in base-pair level models such as the recent work 
by Chou et al.W^ This is presumably due to the 
more complicated structure of the A-form helix 
in dsRNA as opposed to the B-helix in dsDNA. 
However, by comparing to experimental data, 
we have shown that a physical description of the 
properties of dsRNA under torsion and tension 
is still possible. The experimentally observed 
decrease in end-to-end distance with increased 
twist (i.e. positive twist-stretch coupling) of RNA 
is captured well by oxRNA. By contrast, our 
coarse-grained model of DNA does not reproduce 
the anomalous (ne gative ) twist-stretch coupling 
observed in dsDNApZMai note that the model 
of Ref. Uni has reported negative twist-stretch 
coupling for both dsDNA and dsRNA. This 
suggests that, although both positive and negative 
twist-stretch coupling can be represented in the 
framework of coarse-grained models, capturing 
the differential behaviour in both molecules may 
be beyond the scope of present coarse-grained 
descriptions. 


We have investigated the mechanical response 
of dsRNA to twist and stretching force in a 
coarse-grained computational model. To our 
knowledge, this is the first full determination of 
the buckling behaviour of dsRNA in a model 
at single-nucleotide resolution that consistently 
incorporates the salt-dependent thermodynamics 
of double strand denaturation. Reproducing the 
persistence lengths in a quantitatively accurate 
fashion has proven more challenging for dsRNA 
than for dsDNA in the framework of coarse¬ 


Our model is unable to capture the disappear¬ 
ance of the positively supercoiled plectonemic state 
at higher stretching forces. Given the simplified 
nature of the oxRNA model, it is perhaps not too 
surprising that we are unable to capture this “P- 
RNA” state, however we note that the structure 
and physical origins of this state are not yet fully 
understood. 

Similar to our simulations of DNA, we observe 
plectonemes with denaturation bubbles at the tips 

















of their end-loops for negative supercoiling and 
intermediate stretching forces of approximately 
2pN. This coupling of denaturation and writhing 
occurs because the highly bent tip of a plectoneme 
is a particularly favourable location for the nucle- 
ation of a bubble; similarly a bubble is a favourable 
site at which to initiate writhing. In contrast to ds- 
DNA, no end-loop denaturations occurred for pos¬ 
itive supercoiling up to stretching forces of 6pN, 
presumably due to the stronger binding between 
Watson-Crick base pairs in dsRNA. When a plec¬ 
toneme with a tip bubble is present, we predict it 
to be preferentially localised in weak parts of the 
strand sequence, by a mechanism analogous to the 
one described for dsDNA.I^ 

Summing up, we have presented a comprehen¬ 
sive study of dsRNA under torsional and exten- 
sional stress. While reproducing the detailed be¬ 


haviour of the molecule remains a challenge for 
coarse-grained modelling, our findings are in good 
agreement with experimental results and provide 
the basis for capturing the behaviour of more com¬ 
plex RNA structures. 
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Supplementary Material 

S-l. EXTENSION OF THE OXRNA MODEL TO INCLUDE SALT DEPENDENCE 

Following the incorporation of salt-dependent interactions in the oxDNA model of DNA^^, we present 
here a similar extension of the oxRNA model of Ref. m to include salt dependence. We parameterise the 
new interaction in the oxRNA model to reproduce the melting temperatures of RNA duplexes at different 
monovalent (Na+) salt concentrations. The details of the fitting procedure used can be found in Ref. [38l 
The additional term introduced into the oxRNA potential to capture salt effects is of a modified 
Debye-Hiickel form 


where 


f4lectrostatic 




if ^smooth ^ ^ ) 

if ^cut > r ^ ^smooth; 

otherwise. 


,, ^ b-b rr. r\ ( 9 effe)^ exp (-r'" '^/X-DuiT,!)) 


and 


Adh(T, I) 


eperksT 
2NAe^I ■ 


Kimooth is given by 


bsniooth — b {t ^cut) 


(SI) 


(S2) 


(S3) 


(S4) 


with b and rcut chosen so that ^electrostatic is smooth and differentiable. This truncation of Vdh at finite 
distance rcut allows for much faster calculation of forces and pairwise energies between particles. We 
set rsmooth = 3Adh, the same as for the oxDNA2 modelp^ where only negligible differences in oligomer 
melting temperatures were found when using even larger rguiooth- fa the equations above, I is the molar 
salt concentration, e is the electron charge, is the Boltzmann constant, Na is Avogadro’s number, T 
is the temperature, ep is the vacuum permittivity and Cr is the relative permittivity of water (which we 
set to 80). The distance between the interacting sites, which are placed on the backbone sites of the rigid 
bodies representing the nucleotides in oxRNA, is denoted as r^~^. 

In Debye-Hiickel theory, qes is 1. Here, we used the fitting procedure of Ref. [3H|to find the optimal 
value of Qes for the coarse-grained model by fitting it to the melting temperatures of 5, 6, 7, 8, 10 
and 12 mers at salt concentrations ranging from 0.1 M to 0.5 M. The fitting was performed using the 
average-base oxRNA model, to which 14iectrostatic had been added. To obtain the melting temperatures 
of the RNA duplexes to which we fitted the model, we use the melting temperatures as predicted by the 
nearest-neighbour model by Turner et where the respective free-energy contribution of each base 

pair to the duplex stability have been averaged over all possible combinations of Watson-Crick base-pair 
step^^. The nearest-neighbour model was derived for 1 M salt. To obtain the melting temperatures for 
lower salt concentration, we correct the free-energy stability of a duplex by adding an extra destabilizing 
term to the duplex entropy taken from Ref.^ 

AS{N,I) =0Mm\og{I) calmoP^K-i (S5) 


where N is the number of phosphates and / is the molar salt concentration. A duplex can have phosphates 
present at both 3’ and 5’ ends of each strand, but can also have the phosphates cut at one of the ends of 
each strand. As our coarse-grained model does not include an explicit representation of the phosphate 
group, we chose the magnitude of the charges placed on the nucleotides at both the 3' and 5' ends of the 
strand to be (?eff/2. This choice leads to a total charge on the RNA duplex that will be the same as if the 
phosphate charges were cut at one of the ends. Thus, it should be kept in mind that the oxRNA model 
cannot reproduce subtleties caused by having the phosphates cut off one or both ends. 
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Salt concentration [M] 

Motif 

0.1 

0.3 

0.5 

6-mer 

8-mer 

10-mer 

26.3(1.3) / 25.9 
44.9(4.0) / 46.6 
54.6(2.0) / 58.5 

31.2(2.0) / 29.7 
52.9(1.3) / 50.7 
65.4(2.1) / 62.8 

33.6(0.9) / 31.4 
54.2(0.4) / 52.6 
65.4(0.5) / 64.8 


TABLE S-I. The melting temperatures of RNA duplexes at different salt concentrations for the average-base 
parametrization of oxRNA with the new salt-dependent term included (Tm) compared to the melting temperatures 
of the averaged nearest-neighbour model with the salt correction as introduced in Eq. S5 (Tm(NN‘*''®)). The 
individual cells in the table are in the form Tm (error) / where the error was calculated as the 

standard deviation of the melting temperatures estimated from 5 different independent simulations. The melting 
temperatures Tm were estimated from VMMC simulations and are for a strand concentration of 4.2 x 10“® M. 


The correction to the entropy contribution for the nearest-neighbour model in Eq. |S5| is based 
on the hairping unzipping experiments in Ref. 1371 where the stability of a hairpin was obtained 
for varying salt concentrations and temperatures. The average destabilization free-energy was ob¬ 
served to be —0.054Af log (I) kcal/mol, which is similar to that for DNA duplexes at 37°C, which is 
^^37 = —0.057A^log (!) according to Ref. HQI In the nearest-neighbour model for DNA melting in 
Ref. [501 this destabilisation is taken to be only of entropic origin. We hence interpreted the destabi¬ 
lization derived from the RNA hairpin unzipping experiments also as contributing to the entropy in the 
nearest-neighbor model for RNA thermodynamics. If, however, a more detailed study of RNA duplex or 
hairpin thermodynamics at varying salt concentrations becomes available, we might need to revisit our 
parametrization and fit it to more accurate estimations of melting temperatures at varying salts. 

We obtained equal to 1.26 from the fitting procedure. We note the resulting geff is larger than 
1, but given the complexity of potential salt effects, and the simplicity of our mean-field Debye-Hiickel 
representation, not too much can be read into these numerical values. 

To test the fitted value of geff, we studied the melting temperatures of several RNA duplexes at 
varying salt concentrations with virtual-move Monte Carlo simulations (VMMC), using the variant from 
the Appendix of Ref. |5T] Each simulation was run for at least 3 x 10^^ steps. The results are shown in 
Table [S^ for the average-base oxRNA model with the new salt dependence included. 


S-ll. BOUNDARY CONDITIONS 


In order to keep the superhelical density in a dsRNA double strand constant during a simulation, 5 base 
pairs were added to the 600 bp-system at each end, and constrained in stiff, two-dimensional harmonic 
traps. These traps only exert forces in the plane perpendicular to the setup axis of the double strand, thus 
not causing any linear elongation of the system. Analogous constraining boundary conditions have been 
successfully used before in simulations of cruciform extrusiorPSl and dsDNA plectoneme structure^^. A 
schematic overview of the boundary conditions applied is shown in Fig. jSlj 

The two-dimensional harmonic traps used to keep the superhelical density of the system constant are 
implemented by a potential of the form 


1 ■ ' 

Vtrap(rn; ^n,o) = ^ ^ ' ^trap(^n ~ ^n,o) > (^6) 

i-1 

where r„ = (r^, r^, r^) is the centre-of-mass position of the n-th trapped nucleotide and the corresponding 
trap position is r„ g = (I'm gj*"^ gi^n o)> chosen initially such as to fix a given twist angle of the strand. 
We found that choosing = 58.7N/m and kf^'^ = 0 kept the superhelical density fixed by 

preventing rotations of the 5-bp handles at the double strand ends, while not hindering strand extension 
along the setup axis X 3 . 

The RNA duplexes studied in this work have finite length, which means that more distant parts of the 
system can pass around the strand ends. Such a process would modify the superhelical density a of the 
system. Therefore, such movements of the system are prevented in our simulations by repulsion planes 
oriented perpendicular to the setup axis X 3 which co-move with the first boundary nucleotide of the two 
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(a) 


(b) 



FIG. SI. Schematic depiction of the boundary conditions used, illustrated for the last 2 bp at each end of the 
dsRNA system: (a) View along the double strand axis. 5 nucleotides on each strand end are constrained by 
2-dimensional harmonic traps, which fix the boundary nucleotides to positions r„,o in planes perpendicular to the 
strand axis (green), (b) View perpendicular to the double strand. Due to the 2-dimensional traps, nucleotides 
are unconstrained only in the strand-axis direction. A repulsion plane perpendicular to the strand axis is tagged 
to the last base pair. Movement of nucleotides into the area below the end base pair (shaded grey) is therefore 
excluded. In order to allow unconstrained strand extensibility, the repulsion plane does not act on the first two 
base pairs along the strand. 



single RNA strands in the system. Repulsion planes generate a potential 

Vpw(r; R) = ((r - R) • 6)" 9{- (r - R) • 6), (S7) 

where r is the centre-of-mass position of an affected particle, R and 6 are anchor point and orientation 
of the plane, and 6 is the Heaviside step function. We choose 6 = X 3 and 6 = —X 3 for the lower and 
upper repulsion planes respectively, and set R equal to the instantaneous positions of the first and last 
double strand boundary base pair. To avoid restricting free strand extensibility in the X 3 direction, the 
repulsion planes are set up to not interact with the next-to-last boundary base pairs at both strand ends. 
In all simulations, we chose parameters = 29.3pN/nm, which prevented the duplex from passing 

around its ends during all simulation runs. 


S-l 


DETERMINING EXTENSIONAL PROPERTIES OF DSRNA 


The full “hat curves” for all stretching forces F = 0.5,1.0,1.5, 2.0, 2.5,3.0 and 6.0 pN and superhelical 
densities —0.10 < a < -1-0.10 are shown in Fig. |S2[ In order to measure the decrease of end-to-end 
extension as a function of added superhelical density in the post-buckling regime, we fitted linear functions 
to the overtwisted branch of the hatcurves (see Fig. S2). For low stretching forces, the pos tbuckling slope 
of the hat curves decreases at high levels of supercoiling. As has been noted beforel^EIMl^ tbis finite-size 
effect is due to the interactions of the double strand with the system boundaries. In order to obtain 
the generic behaviour of the system, we attempted to restrict the fitting to a range in a in which the 
postbuckling curve exhibits no non-linearities (see Fig. S2). There is some ambiguity in choosing the 
range of the linear fits. We therefore performed two separate fits where we shifted the fitting domain by 
one point towards the buckling transition, as shown in Fig. S2 The values shown in Fig. 3(c) of the main 


paper refer to the mean and standard deviation of the two values obtained in this way. 

The slopes thus determined can then be directly compared to experimental results, as shown in Fig. 3(c) 
of the main text. 

The measured values of the post-buckling slopes can furthermore be used to determine the twist stiffness 
of the plectonemic state P by fitting to a relation obtained by MarkcP^. Following Refs. [22l and ITBI this 
relation is: 


dAL 


Po 


1-5 


keT 

AoF 


16 


r^y 

AqF ) \c\l 1-p/c J 


/ 2pg / 1 _ 

V 1-P/c VP 


dALk 


(S 8 ) 
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FIG. S2. Mean strand end-to-end extensions for all values of parameters F and cr studied in this work. Postbuck- 
ling slopes were determined by fitting linear functions to the overtwisted branches of the hat curves. As the slopes 
obtained somewhat depend on the fit range chosen, we performed two different hts over slightly different ranges 
of (T, as shown in (a) and (b). Using a thermodynamic model due to MarkcP^, the stiffness of the plectonemic 
phase can be determined from this post-buckling slope. 


where Aq and Cq are the bending and twist persistence lengths, po is the equilibrium helical pitch, 
00 the equilibrium twist angle and F the stretching force. Furthermore, g = F — yjFksT/A q, wh ile 
p = HeTPOq and c = ksTCOl are proportional to P and Cq respectively. The result of fitting Eq. S8 to 
the postbuckling slopes determined from simulations is shown in Fig. 3(c) of the main text. 


S-IV. DETECTION OF DOUBLE STRAND MELTING AND PLECTONEME POSITION 

As the value of Vhb paired nucleotides assumes continuous values, it is necessary to dehne a cutoff 
criterion to dete rmine whether a given pair of nucleotides is base-paired or not. Following the approach 
taken previous!}®!!!]^ counted a base-pair as formed when the interaction energy from hydrogen 
bonding between two nucleotides was below —4.13 x 10“^^ J, corresponding to approximately 15% of the 
typical energy of a fully formed hydrogen bond. 

In order to assign a position variable to a given plectoneme structure, we used the plectoneme detection 
algorithm described in detail in Ref. HSl The individual steps of plectoneme detection ar^^ 

• Start from a double strand end, loop over all base pair centre points 

— If any part of the remaining double strand that is more than bp away along the contour of 
the duplex has a distance dim < record the index of the current base pair as the beginning 
of a plectoneme, if the beginning of a plectoneme has not yet been detected before. 

— If diin > d°jjj for all base pair centres of the remaining double strand and a plectoneme beginning 
has been detected before, record the current base pair index as the end of a plectonemic region 
and continue searching for further plectonemes from the next base pair centre 

• The plectoneme position is the mean between the base pair indices of the beginning and end of a 
plectonemic region 

• The plectoneme size is the difference between the base pair indices of the beginning and end of a 
plectonemic region 

The systems studied in the present work are simulated at a monovalent ionic strength of 100 mM, 
which is significantly lower than the 500 mM ionic strength considered for the analogous dsDNA system 
in Ref. [281 As a consequence of the increased electrostatic strand repulsion due to lower salt, the diameter 
of the end-loop and plectoneme stem are expected to slightly increase. It was found that the properties of 
these somewhat larger structures is best captured when setting the detector parameters to = 10.1 nm 
and Nc = 50 bp, which are slightly larger than the values used in Ref. |2H1 We note that W represents 
a lower limit on the size of plectoneme structures that can be detected using the detection algorithm 
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outlined above. However, at an ionic strength of 100 mM, typical plectoneme structures are significantly 
larger than 50 bp, and are therefore reliably detected by the algorithm. 

As in our previous work on dsDNA ('Ref.[^5|l. a tip-bubble plectoneme is defined as a plectoneme whose 
midpoint as dehned by the detection algorithm is less than 20 bp away from the centre of a denaturation 
bubble. 
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